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1 .  Introduct i on 

^The  objective  of  the  research  undertaken  was  to  provide  a  theory  of 
consolidation  of  fine-grained  soils  of  sufficient  generality  to  enable 
predictions  to  be  made  in  engineering  situations  where  the  sediment  is  so 
soft  that  allowance  must  be  made  for  changes  in  material  properties 
during  the  progress  of  consolidation,  where  large  strains  and 
displacements  must  occur,  account  being  taken  of  pore  water  flow  in  two 
or  three  dimensions. 

The  study  proceeded  in  two  Stages.  The  first  Stage  involved  a 
critical  examination  of  presently  available  theories  of  consol idat ion 
which  were  modified  and  extended  to  meet  the  objective  of  the  research. 
The  aim  in  Stage  I  was  to  develop  mathematical  equations  governing  the 
consolidation  process  based  on  physical  assumptions  which  accord  with  the 
known  behavior  of  soft  sediments. 

In  Stage  II  numerical  procedures  and  associated  algorithms  were 
considered  based  on  the  analysis  which  would  permit  prediction  of  the 
distribution  and  time  variation  within  the  medium  of  quantities  of 
engineering  interest,  namely:  the  pore  water  pressure,  the  effective 
stresses  and  the  void  ratio. 


2.  Theories  of  Consolidation:  A  Brief  Survey 


The  mechanism  of  the  consolidation  process  in  fine-grained  soils  and 
its  expression  in  analytical  form  were  given  first  by  Terzaghi  in  1923 
(*).  The  theory  was  restricted  to  pore  water  flow  in  one  (Cartesian) 
dimension  but  as  was  recently  pointed  out  did  take  some  account  of  large 
strains  (  ).  It  was  aimed  at  providing  a  means  for  predicting  the 
progress  of  settlement  in  uniformly  loaded  clay  layers.  Since  that  time 
the  theory  of  one-dimensional  consolidation  has  been  greatly  extended  to 
take  account  of  non-linear  behavior,  including  the  variation  during 
consolidation  of  the  coefficients  of  permeability  and  compressibi l ity, 
and  other  effects  associated  with  large  strain  and  deformation  (see,  for 
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example,  3,  ** ,  5)  and  this  development  was  made  possible  by  the  intrinsic 
simplicity  of  one-dimensional  deformation. 

The  problems  of  consolidation  which  face  the  civil  engineer  are 
rarely  one-dimensional  and  a  need  soon  developed  for  a  more  general 
formulation  of  the  theory  which  would,  for  example,  allow  the  settlement 
of  structures  to  be  predicted  or  sand  drain  installations  to  be 
rationally  designed:  both  cases  where  soil  deformation  and  pore  water 
flow  occur  in  three  dimensions.  Revere  difficulties  were  soon 
encountered  in  this  quest  and  it  was  found  that  only  by  introducing 
radical  simplifying  assumptions  could  any  progress  be  made.  This  phase 
of  development  is  associated  with  Terzaghi  and  Rendulic  -  both  engineers 
-  and  with  the  mathematician  Maurice  Biot.  The  nature  of  their 
contributions  reflect  their  difference  in  outlook  and  objective.  Biot 
developed  a  rigorous  and  self-consistent  theory  (6)  based  on  the 
assumption  that  the  soil  skeleton  behaves  as  a  porous  perfectly  elastic 
medium,  although  this  is  only  a  crude  approximation  to  the  known  behavior 
of  real  soils.  Terzaghi,  on  the  other  hand,  tried  to  avoid  commitment  to 
particular  constitutive  relations  between  components  of  strain  and 
effective  stress  and  he  succeeded,  but  only  at  the  cost  of  implying  that 
during  consolidation  (under  constant  loading)  all  components  of  total 
stress  remain  constant  everywhere  (  ,  ).  In  reviewing  these  theories  in 
1943  (9)  Terzaghi  remarks  "...the  existing  theoretical  methods  for 
dealing  with  two-  and  three-dimensional  problems  of  the  consolidation  of 
clay  under  load  are  not  yet  ready  for  practical  application".  Owing  to 
the  complexity  of  the  problem  progress  was  slow,  but  by  1969  the 
consequences  of  these  two  differing  approaches  had  been  examined  in  a 
sufficient  number  of  particular  cases  for  some  general  conclusions  to  be 
drawn  by  Schiffman  and  co-workers  (*®). 

In  all  these  studies  the  assumption  was  invariably  made,  either 
tacitly  or  explicitly,  that  during  consolidation  the  effective  stress- 
strain  relations  are  linear,  the  coefficient  of  permeability  of  the  soil 
remains  constant  and  that  the  strains  and  deformations  of  the  soil  are 
small.  However,  these  factors  can  have  an  important  influence  on  the 


progress  of  consolidation,  in  particular  upon  the  magnitude  and  changing 
pattern  of  pore  water  pressure  distribution.  More  recently  engineers 
have  sought  to  remove  those  restrictions  from  the  theory,  but  owing  to 
the  complexity  of  the  general  case  have  confined  their  attention  to 
problems  of  one-dimensional  compression  and  flow  (see,  for  example, 

18).  Contemporaneously,  Biot's  work  has  been  extended  to  large  strains 
and  deformations  (  ,  ,  ,  ),  but  the  "soil"  considered  remains  highly 

ideal i zed .* 

3.  The  Present  Study  . 

The  typical  problem  with  which  we  are  concerned  here  is  the 
consolidation  under  its  own  weight  of  a  mound  of  soft  saturated  clay 
which  has  been  placed  beneath  water.  It  has  been  assumed  that  the 
strains  and  deformations  that  may  develop  during  consolidation  are  large: 
that  is,  account  must  be  taken  in  the  theory  of  the  variation  of 
compressibility  and  permeability  during  consolidation. 

Attention  has  been  restricted  to  two-dimensional  pore-water  flow  and 
one-dimensional  soil  deformation.  The  details  of  the  analyses  are  given 
in  Appendixes  A  (Eulerian  description)  and  B  (Lagrangian  description). 
Appendix  C  formulates  the  governing  relationships  for  an  axially 
symmetric  geometry  using  a  Lagrangian  description.  The  most  promising 
approach  seems  at  present  to  lie  with  a  Lagrangian  description  allied 
with  the  void  ratio  (e)  as  the  dependent  variable.  The  governing 
equations  ( B 1 9 ,  B26)  which  result  are  non-linear  and  must  be  solved 
numerically  using  an  iterative  procedure. 

A  numerical  study  of  specific  cases  may  reveal  that  strains  and 
deformations  are  not  unduly  large,  in  which  case  alternative 
approximations  to  those  adopted  may  prove  to  be  more  appropriate.  This 
is  likely  to  be  helpful  if  the  more  general  problem  of  successive 

*We  remark  that  a  theory  of  large  strain  consolidation  which 
regards  the  coefficient  of  permeability  as  a  soil  constant  is  a 
contribution  exclusively  to  Applied  Mechanics. 


accret  ion  of  material  on  an  existing  mound  is  to  be  treated  economically, 
although  the  present  theory  can,  in  fact,  be  used  to  cover  this  case. 

4 .  Numerical  Analysis:  One-Dimensional  Theory 

Several  numerical  analysis  techniques  can  and  have  been  applied  to 
the  equation  governing  one-dimensional  consolidation.  These  include  an 
implicit  finite  difference  technique  (2^),  and  explicit  finite  difference 
technique  (2I*)*  and  the  method  of  lines  (2^).  The  finite  difference 
procedures  approximate  the  spatial  derivatives  by  centered  differences 
and  the  time  derivative  by  a  forward  difference.  The  method  of  lines  is 
a  classical  method  which  reduces  a  partial  differential  equation  to  a 
system  of  ordinary  differentia!  equations  by  discretizing  all  but  one  of 
the  independent  variables.  Given  a  partial  differential  equation  of  the 
form  of  eauation  (A22)  and  using  a  centered  difference  approximation  for 
the  spatial  coordinate,  one  obtains  a  system  of  ordinary  differential 
equations  of  the  form 

8u  i 

yjT"  =  f(t,u1,u2,  ...  u^):  i  =  1,2, ...,m  (1) 

where  thpre  are  m  spatial  mesh  points. 

The  system  of  ordinary  differential  equations  is  solved  by  an 
implicit  technique.  An  Adams-Bashforth  predictor  is  used  in  coniunction 
with  an  Adams-Moul ton  corrector.  If  the  system  of  equations  (1)  is 
stiff,  the  predictor-corrector  method  is  coupled  with  Newton's  method  to 
solve  the  nonlinear  equations  by  iteration  (26). 

A  number  of  software  packages  exist  for  the  purpose  of  solving  one¬ 
dimensional  nonlinear  finite  strain  consolidation  problems. 

•  A  suite  of  packages  have  been  developed  by  Bromwell  Engineering, 
Inc.  (now  Bromwell  and  Carrier,  Inc.)  of  Lakeland,  Florida, 
these  packages  are  based  on  the  governing  equation  expressed  in 
terms  of  the  excess  pore  water  pressure.  They  use  an  implicit 

*We  prefer  the  explicit  technique  since,  unlike  the  implicit 
procedure,  a  stability  criterion  exists. 


finite  difference  solution  technique.  These  programs  assume  that 
the  void  ratio-permeability  and  void  ratio-effective  stress 
relationships  are  governed  by  power  laws.  They  have  been  used  in 
practice  in  the  design  and  analysis  of  dredged  fills  and  mine 
waste  planning  (23,2^,28). 

A  series  of  programs  have  been  developed  at  the  University  of 
Colorado  and  at  the  Waterways  Experiment  Station.  These  programs 
have  been  based  on  the  assumption  that  the  void  ratio-effective 
stress  and  void  ratio-permeability  relationships  are  exponential 
functions.  These  programs  are  based  on  a  governing  relationship 
where  the  void  ratio  is  the  dependent  variable.  An  explicit 
finite  difference  procedure  has  been  employed  to  solve  the 
governing  equation.  The  programs  developed  at  the  University  of 
Colorado  were  developed  primarily  for  research  purposes  (  ). 

The  second  generation  programs  developed  at  the  Waterways 
Experiment  Station  were  designed  primarily  to  assist  in  the 
planning  of  dredge  fill  operations  (29,30). 

The  University  of  Colorado  has  developed  a  series  of  programs 
which  are  unrestricted  with  respect  to  the  form  of  the  void 
ratio-affective  stress  and  void  rat io-perraeabi 1 ity  relationships. 
These  programs  are  based  on  the  void  ratio  as  the  dependent 
variable  in  the  governing  equation.  They  use  the  method  of  lines 
as  the  numerical  procedure  of  choice.  They  have  been  used  for 
research  purposes  and  to  analyze  mine  waste  and  marine 
geot echoo logy  problems  (»>»>>>»'* 

A  program  developed  for  Ardaman  and  Associates,  Orlando,  Florida 
has  been  reported  (38>.  Unfortunately,  at  the  time  of  this 
writing,  there  are  insufficient  details  available  concerning  the 
algorithms,  numerical  method  and  use  of  this  software  item. 


Numerical  Analysis:  Two-Dimensional  Theor 


There  are  four  numerical  procedures  which  can  be  applied  to  the 
governing  equation  for  two-dimensional  nonlinear  finite  strain 
consolidation  developed  in  the  Appendices. 


First,  an  explicit  finite  difference  procedure  can  be  developed.  In 
this  procedure  the  spatial  derivatives  can  be  expanded  in  terms  of 
centered  differences,  at  time  t.  The  time  derivative  can  be  expanded  as 
a  forward  difference.  This  will  result  in  an  explicit  difference 
equation.  The  time  marching  procedure  should  be  investigated  for 
stability.  However  stability  will  be  assured  if  each  mesh  size  Aa  and 
Ab  and  the  time  increment  At  meets  the  stability  criterion  for  the 
two-dimens ional  problem  at  each  time  step. 


The  second  procedure  is  implicit.  Here  the  spatial  derivatives  are 
expanded  in  terras  of  centered  differences  at  time  (t  +  At).  This  pro¬ 
cedure  will  establish  a  set  of  simultaneous  equations  for  which  an 
alternating-direction  implicit  procedure  may  be  feasible.  However,  the 
stability  criterion  for  implicit  formulations  is  not  established. 

The  third  suggested  procedure  is  a  variation  of  the  method  of  lines 
described  previously.  Here  a  system  of  ordinary  differential  equations 
of  the  form  (see  Appendix  D) 


3e .  . 

=  f(t,e.^); 


1,2,...  ,m 

1  >2 . n 


(2) 


is  obtained.  The  one-dimensional  results  in  ra  ordinary  differential 
equations.  This  two-dimensional  problem  results  in  m^n  ordinary 
differential  equations.  These  equations  can  be  solved  by  the  same 
techniques  as  discussed  for  the  one-dimensional  formulation.  This  method 
is  always  stable. 


A  fourth  possible  method  would  apply  a  finite  element  discretization 
to  the  spatial  variables  along  with  an  explicit  time  marching  scheme. 

This  approach  has  the  disadvantage  that  a  variational  principle  for  this 
class  of  problems  must  be  established  a  priori. 


It  is  noted  that  the  governing  (B19)  equation  contains  an  integral 
term  from  equation  (B25).  This  can  be  discretized  by  means  of  Simpson's 


rule . 


The  suggested  solution  algorithm  is  as  follows: 

•  Given  the  boundary  vilue  problem  in  terms  of  the  void  ratio  a 
solution  is  developed  for  e(a,b,t). 

•  The  vertical  effective  stress  a'(a,b,t)  is  determined  from  the 
void  rat io-ef fective  stress  relationship  which  is  known  a  priori. 

•  The  vertical  total  stress  o(a,b,t)  is  determined  from  the  bulk 
equilibrium  relationship 


[npf  +  (l-n)p  ]  [j~]  =  0 


•  The  static  pore  water  pressure  uo(a,b,t)  is  determined  by 
integrating  the  fluid  equilibrium  relationship, 

3u  . 

alT  +  T+e~  pf  =  0  < 

o 

(We  note  that  this  assumes  a  constant  level  of  height  of  water 
above  the  base  of  the  clay  layer.  The  consideration  of  a 
variable  free  water  surface  is  a  trivial  extension  to  this 
algorithm. ) 

•  Fach  a  station  at  each  time  of  interest  is  converted  to  the 
appropriate  5  station  by 


t  (.  l  .  fS  l+e(a' ,b,t)  .  , 
5  (a,b,t)  -  /  1+e(a. >b>0)  da 


The  development  of  a  detailed  step-by-step  algorithm  and  the  code  to 
implement  the  algorithm  are  beyond  the  scope  of  the  current  contract. 
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APPENDIX  A 

In  this  Appendix  we  consider  the  case  of  pore-water  flow  and  soil 
skeleton  deformation  restricted  to  two  dimensions  (the  x,z-plane)  and  we 
seek  to  derive  the  equations  governing  the  motion  of  these  phases  by 
adopting  the  so-called  Eulerian  scheme  of  description.  We  therefore 
consider  the  soil  grains  and  water  which  at  any  time  (t)  reside  within 
the  rectangular  element  of  sides  (dx,  dz)  located  at  the  point  (x,z). 

We  denote  by  (w^,  w  )  the  velocity  components  of  the  solid  phase  and 
by  (v^,  v^)  those  of  the  pore  water.  These  components  will  be  functions 
of  the  three  independent  variables  (x,z,t). 


I .  Equations  of  Continuity 

The  volume  porosity  (volume  of  void  space  V  per  unit  of  bulk 
volume  V)  we  denote  by  n,  so  that 


n  =  V  /V  , 
v 


while  the  area  porosity  (area  of  void  space  Av  per  unit  of  bulk  area 
A)  we  denote  by  n  .  If  the  number  of  grains  occupying  the  rectangle 
is  in  some  sense  large,  it  can  be  shown  that 


n  =  n 
a 


The  rate  of  increase  of  weight  of  solids  within  (<5x,  <5z)  must  equal  the 
net  rate  of  flow  of  weight  across  the  four  faces  of  the  element  and  this 
leads  us  to  the  equation 

[  ( 1-n) o  ]  +  [w  (l-n)p  ]  +  4—  [w  (l-n)p  ]  =  0  (All 

3t  c  s J  9x  1  x  sJ  3z  1  z  sJ 


where  p  is  the  unit  weight  of  the  solids, 
s 

Similarly,  for  the  pore-water 


V-j 


where  is  the  unit  weight  of  the  fluid.  These  are  the  equations  of 

cont inui ty . 

Since,  in  the  range  of  stress  with  which  we  are  concerned,  the  unit 
weights  Pg,  p^  can  be  regarded  as  constants,  it  follows  from  (Al)  and 
(A2)  by  addition,  that 

I—  f nv  +  (l-n)w  ]  +  - —  [nv  +  (l-n)w  ]  =  0  .  (A3 

dx  1  X  XJ  o Z  L  Z 


xJ  dz  L  z 


2 .  Flow  Rule 

We  shall  assume  that  the  pore-water  moves  through  the  soil  skeleton 
in  accordance  with  Darcy's  law,  due  account  being  taken  of  the  fact  that 
it  is  the  relative  velocity  between  the  phases  which  induces  a  drag  on 
the  soil  skeleton.  Accordingly,  we  write 


,  ,  _  x  3u 

n(v  -  W  )  - - -r— 

x  x  Pj  3x 


,  ,  z  3u 

n(v  -  w  )  =  -  - — 
z  z  Pj.  3z 


where  (k^,  k  )  are  the  coefficients  of  permeability  appropriate  to  the 
(x,z)  directions  and  u  is  an  excess  pore  water  pressure  which  is 
defined  in  terms  of  the  pore-water  pressure  (p)  by 


u  =  p  +  P^z  +  const. 


where  the  positive  direction  of  z  is  against  gravity. 


3.  The  Governing  F.quations 


If  we  set  (A& )  ,  (A3)  in  (A3)  we  find  that 


3 


and  it  can  be  seen  that  if  the  soil  skeleton  is  stationary  (w^  =  0, 

w  =0),  then  a  familiar  equation  governing  the  excess  pore-water 
z 

pressure  is  obtained,  namely 

3  r.  3ui  .  3  r._  3'»i  _  n 


[k  ~]  *  f-  [k  =  o  . 

3x  L  x  3xJ  3z  1  z  3zJ 


It  is  worth  noting  that  iust  after  dumping  a  mound,  w^  =  0  and  w^  =  0 
almost  everywhere  and  so  (A8)  holds  at  this  instant:  this  eauation  again 
holds  when  consolidation  is  complete. 


Now,  from  (A1 ) 
3w 


1  f  3n  3n  3nl 

rral  3t  +  Wx  3x  +  wz  3zJ 


3x  3z  ( 1-n)  {_  3t  x  3x  z  dz J 

and  eliminating  the  rate  of  soil  skeleton  dilatation  between  (A7)  and 
(A9)  we  find 


3__ 

-  k 

x 

3u 

B 

k 

z 

Bu 

Bx 

-  pf 

Bx_ 

Bz 

-Pf 

Bz " 

))]  =  0 


(A10) 


where  the  operator 


D  _  3  3  3 

- -  -  4-  w  - 4-  W  —  . 

Dt  ”  3 1  x  3x  z  3z 


In  order  to  proceed  further  the  relations  connecting  strain  increments 
with  the  effective  stresses  and  the  effective  stress  increments  must  be 
assumed.  These  relations,  together  with  the  equations  of  equilibrium, 
would  allow  a  complete  theory  to  be  developed.  However,  in  view  of  the 
wide  varietv  of  soil  types  and  their  complicated  response  to  load,  a 
completely  general  approach  is  at  present  not  feasible. 


We  shall  therefore  proceed  in  the  spirit  of  Terzaghi  and  Rendulic's 
early  work  and  seek  plausible  simplifications  which  will  permit  solutions 
to  be  obtained  to  problems  which  will  be  sufficiently  exact  for 
engineering  purposes. 


4 .  One-Dimensional  Compression 

We  shall  assume  here  that  the  development  of  a  mound  takes  place  in 
a  number  of  stages  each  consisting  of  the  sudden  dumping  of  soil  followed 
by  a  period  of  consolidation  during  which  no  further  loading  occurs. 

During  the  consolidation  it  will  be  assumed  that: 

(a)  Lateral  displacement  of  the  soil  can  be  ignored. 

The  mound  will  therefore  consolidate  with  the  development  of 
vertical  strains  only  but  the  pore  water  flow  we  shall  not  assume  to  be 
constrained  in  any  way. 

The  porosity  (or  void  ratio)  at  any  point  in  the  clay  will  depend  on 
the  initial  porosity  and  state  of  effective  stress  there  and  also  upon 
the  subsequent  development  of  the  components  of  effective  stress  as 
consolidation  proceeds.  The  detailed  effective  stress-strain  behavior  of 
soils  is  very  complicated  and  to  obtain  an  engineering  solution  to  the 
problem  we  introduce  the  following  simplifying  assumptions. 

(b)  The  void  ratio  (e)  of  a  soft  saturated  clay  depends  to  a  good 

approximation  only  upon  the  major  principal  effective  stress  cx ^ ' 
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(the  so-called  "American  Hypothesis";  see,  for  example,  ref.  “  ). 

(c)  The  major  principal  stress  is  (almost  everywhere)  equal  to  the 

vertical  stress  a  (see  Ref.  16)  in  a  mound-like  structure.  It 
zz 

follows  from  these  two  assumptions  that 

e  =  e( a  ' )  =  e( a  ').  (All) 

1  zz 

For  simplicity  of  exposition  we  restrict  our  remarks  to  the  case  where 
k^  =  k^  =  k(const);  the  more  general  problem  can  if  required  readily  be 
treated  later  in  detail.  Equation  (A10)  can  be  written  in  the  form: 

-15- 


where 


Since 


k_  2  =  1  de  _ 

Pf  U  ( 1+e)  dOj 7  Dt 


D  _  3  3_ 

Dt  3t  Wz  3z 


(A12) 


(A13) 


which  is,  apart  from  the  factor  (1+e)  in  place  of  (1+e^),  Terzaghi's 
coefficient  of  one-dimensional  consolidation. 


5.  One-Dimensional  Consolidation 

It  is  worth  noting  that  in  cases  of  strictly  one-dimensional 
compression  and  pore  water  flow,  the  above  assumptions  (a),  (b)  and  (c) 
are  satisfied  and  equation  (A14) 


c 

v 


,2  _  Da 

3  p  Dp  zz 

-2  Dt  Dt 

3z 


(A15) 


is  exact. 

This  equation  governing  the  pore  water  pressure  has  a  structure  very 
similar  to  that  encountered  in  Terzaghi's  theory  and  indeed  his  equation 
can  be  recovered  by  replacing  the  differential  operator  D/Dt  by  3/3t. 
From  (A13)  this  can  be  seen  to  be  tantamount  to  ignoring,  for  example, 
w  3p/3*  compared  with  3p/3t. 

There  is  no  a  priori  reason  to  suppose  that  this  can  be  justified. 
Some  indication  of  the  error  involved  can  be  found  from  (A3)  which 
reduces  to 


-16- 


when  there  exists  a  plane  on  which  the  condition  =  w^  =  0  persists. 
(This  is  so  in  the  oedometer  test:  the  mid-plane  with  two-way  drainage, 
or  Che  base  if  there  is  no  flow  from  it.)  It  then  follows  from  (A16)  and 
(.V>)  that 


k  3u 
Wz  Pj  3z 

and  so,  for  example: 

Du  _  £u.  +  /  3u 

Dt  3t  p ^  \3z / 


( Al  7 ) 


What  then  is  the  equation  governing  the  pore-water  pressure  which  is 
exact  within  the  context  of  the  assumptions  made  in  the  theory,  when  the 
Eulerian  description  is  used? 

We  commence  from  the  form  of  (A15)  which  allows  for  the  variation  of 


permeability,  namely 

1  3u 

f  k  3u ] 

_  2£ 

Da 

ZZ 

(A18) 

m  3z 

V 

[Pf  3ZJ 

Dt 

Dt 

where,  from  (A6): 

Dp  _  3u 
Dt  =  Tt 

+  Js-  |“  [ 

pf  3z  [ 

3u 

3z 

’<]• 

(A19) 

Also, 

Da 

zz 

Dt 

3a  .  , 

zz  k  au 

dt  *  3z 

3a 

zz 

3z 

(A20) 

but  vertical  equilibrium  requires  that 


Setting  (A21)  in  (A20)  and  subtracting  the  resulting  equation  from  (A19), 
it  is  found  that  (A18)  becomes 

2  k(p  -pJ  „  „  3o 


Using  (A6)  a  similar  equation  can  be  derived  for  the  pore-water  pressure 

(p). 

This  equation  is  highly  non-linear  even  in  the  "thin"  layer  case 
(p  =  P^.),  and  its  use  in  problems  of  large  displacement  and  strain  is 
limited  for  the  following  reasons: 

(i)  The  parameters  k  and  m^  are  related  to  the  void  ratio  e,  but 
this  in  turn  is  not  connected  to  u  (or  p)  in  a  straightforward 
way . 

(ii)  The  vertical  total  stress  a  and  its  time  derivative  3o  /3t 

zz  zz 

for  a  fixed  value  of  z  varies  with  time  in  a  way  which  is  not 
given  as  part  of  the  data  of  the  problem  and  must  be  discovered  as 
part  of  the  solution: 

(iii)  The  geometry  of  the  mound  boundary  on  which  a  condition  of  the 

type  u  =  0  persists  is,  again,  not  known  ab^  initio  and  must  be 
found  during  the  solution  process. 

However,  one  distinct  advantage  of  working  in  terms  of  u  is  that  the 
initial  conditions  can  be  specified  in  terms  of  the  initial  state  of 
total  stress  in  the  mound.  This  allows  a  more  complete  and  exact 
description  than  can  be  achieved  when  the  void  ratio  or  porosity  is  used 
as  the  independent  variable  (see  Appendix  B) . 

A  rather  more  promising  governing  equation  emerges  if  the  porosity 
(n)  is  taken  as  the  dependent  variable.  Commencing  from  (A10)  we  can 


write 


so  that  using  (A21) 


3u 

3z 


do' 


=  -  ( 1-n) (p  -p  )  - 

s  r  dn 


3n 

3z 


(A25) 


Setting  this  expression  into  (A23)  we  find,  after  some  algebra,  that 


ll  tcv  &  ■  *  [^7  -  ']  It  Ifc(l-n)2] 


3n 

3z 


(A26) 


which  is  essentially  equation  (23)  of  (5).  It  is  not  open  to  objections 


(i)  and  (ii)  above.  The  difficulty  (iii)  remains  but  it  has  been  claimed 


that  this  can  be  overcome  by  using  suitable  numerical  techniques*. 


*Lee  and  Sills  (Ref.  5)  use  a  numerical  technique  originating  from 
work  by  Crank  and  Gupta  (^1,^2)  which  effectively  updates  the  position  of 
the  moving  boundary  as  the  solution  proceeds.  They  take  the  case  of  a 
"thin"  soil  layer  where  the  approximat ion  Pg  =  P^  can  be  justified  on 

physical  grounds,  and  compare  their  numerical  solution  with  that  given  by 
Gibson,  England  and  Hussey  (Ref.  4).  The  agreement  is  apparently 
excellent,  but  is  not  wholly  convincing  as  the  coefficient  c  in  ref.  4 


is  taken  incorrectly  as 


2  2  2 
c_  =  c  /(1+e)  instead  of  c„  =  c  (1+e  )  /(1+e) 

V  v  K  v 


APPENDIX  B 


In  this  Appendix  we  consider  the  physical  problem  discussed  in 
Appendix  A,  in  particular  the  case  of  two-dimensional  pore-water  flow  and 
one-dimensional  compression  in  a  submarine  clay  mound.  However,  here  we 
use  the  Lagrangxan  scheme  of  description  and  examine  the  motion  of  an 
element  of  soil  which  contains  the  same  solids  throughout  its  history. 

For  reasons  mentioned  in  Appendix  A  we  shall  not  seek  to  derive  an 
equation  governing  the  pore-water  pressure  (or  an  excess  pore-water 
pressure),  but  rather  work  in  terms  of  the  void  ratio  (e)  or  the  porosity 
(n).  The  same  physical  assumptions  used  there  will  again  be  adopted 
here.  The  argument  which  ensues  follows  closely  that  used  in  ref.  **  and 
we  shall  not  repeat  the  details  here. 


1.  Equations  of  Continuity 

The  deformation  of  the  soil  skeleton  is  described  by  following  the 
motion  of  an  element  of  initial  area  (6a  x  6b)  lying  at  t  =  0  at 
(a,b)*  which  during  its  subsequent  motion  always  contains  the  same 
solids.  The  location  of  a  material  point  initially  at  (a,b)  is  given 
subsequently  by 


5  =  5(a,b,t) 

n  =  n(a,b,t ) 

where,  by  definition, 

a  =  5(a,b,0) 
b  =  n(a , b ,0) 


(Bl) 

(B2) 


(  B  3 ) 
(B4) 


The  coordinate  a 


points  upward  against  gravity. 
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(Bl) 

(B2) 


(B3) 

(BA) 


*The  coordinate  a  points  upward  against  gravity. 


Since  we  shall  assume  that  the  motion  of  the  soil  skeleton  is  confined  to 
the  a  direction  it  follows  that  (B2)  may  be  replaced  by 


n  =  b  . 

The  equation  of  continuity  of  solids  takes  the  simple  form 

H  =  =  1+e 

3a  1-n  1+e^ 

where 

nQ  *  n(a ,b,0) 
eR  =  e(a,b,0) 


(B2)bis 


(B5) 


are,  respectively,  the  initial  (t  =  0)  porosity  and  void  ratio  at  (a.b). 


As  the  element  translates  and  deforms  during  its  subsequent  motion, 
pore-water  moves  into  and  out  from  its  four  faces  and  the  equation  of 
continuity  of  this  (incompressible)  phase  is  found  to  be 


(B6) 


2.  Flow  Rule 


Darcy's  Law  takes  the  form* 


,  ,  k  3u 

niv  -w  )  =  -  —  -r— ■ 
a  a  p  35 


,  ,  k  3u 

nlv  -w,  )  =  -  —  -r~ 
b  b  p  3n 


with  this  description,  but  since 


3u  _  3u  35 

Ja  ~  8a  ’ 


We  have  taken  k 

2 

anisotropic  case  is  trivial 


k^  =  k  here,  but  the  extension  to  the 


(B7) 

(B8) 
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and 

it  follows  that 


3u  3u 

■SF  Tn 

w,  *  0 
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n(v  -w  )  -p- 
a  a  3a 
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k_  3u 
Pf  3a 

k_  3u_ 
Pf  3b 
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(Bll) 
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(B13) 


The  excess  pore-water  pressure  (u)  is  related  to  the  pore-water  pressure 


by 

u 

=  p  +  -  const. 

The  Governing  Equations 

Setting  (B12)  and  ( Bl 3) 

in  (B6)  we  find 

3  fl.  3u  /  l*eo\~ 

+  U{&L  \  k  -^1 

Saj^Pj  3a  y 1+e  J 

3b[\1+eojpf  3bJ 

(B14) 


(B15) 


where  (B5)  has  been  used.  Using  now  (B5)  and  (B14),  equation  (BL5)  can 
be  expressed  in  terms  of  the  pore-water  pressure: 


3_ 

3a 


i_  +  J  +  L.\ *2.  (ill  VI  = 

3a  \  1+e  J  k  3b  I  pf  lb  ^1+eJJ  TT^ 


3e 
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(B16) 


Vertical  equilibrium  requires  that 

!-  +  [°Pf  +  (i-n)%]  H  =  °  > 


(B17) 


where  n  is  the  vertical  total  stress 
conjunction  with 


and  this  may  be  used  in 


which  defines  the  vertical  effective  stress  (o'),  to  eliminate  p  from 
(B16). 


The  resulting  equation  governing  the  void  ratio  (e)  ;s 

3 _ 

3  a 


a  _  r  k_  /  l  +  e  \  3o]  =  1  _3e 

3b[pf  ^1+eo/  3t 


3  f  (l+e)2  3e" 

3b  L  F  (1+eJ3  3bJ' 
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d  T  k  1  3e 

de  L  l+e  J  3a 


in  the  general  case,  which  reduces  when  is  a  constant  to 


*fc[ 


de 

TZ 


where  (as  in  Ref,  4 )  the  coefficient  of  finite  consolidation 


2 


k  (,+V  do1 


CF  (l+e)  de 


(B19) 


(B20) 


and 


GS  "  Ps/Pf  ’ 


(B21) 


4.  The  Total  Stress  Term  * 

The  last  term  on  the  Left-hand  side  of  (B19)  is  not  known  a  priori 
and  must  be  found  and  updated  during  the  solution.  In  this  Section  we 
show  how  the  term  3o/3b  can  be  found  at  any  time  from  the  current  void 
ratio  distribution  e(a,b,t). 

Consider  a  vertical  cylinder  of  material,  of  unit  cross-sectional 
area,  extending  from  the  base  of  the  mound  to  the  water  level  above  the 
mound  and  denote  this  height  by  H.  It  is  evident  that  the  weight  of 
material  (solids  and  water)  within  this  cylinder  would  remain  constant  if 
no  water  flowed  across  the  surfaces  (i.e.,  strict  one-dimensional 
consolidation)  or  was  added  (H  =  constant).  Under  these  conditions 
3o/3b  on  the  base  depends  on  b  alone.  However,  on  the  surface  of 


the  mound  a  =  a^Cb)  settlement  increases  the  water  pressure  and  hence 
the  total  vertical  stress  there.  Thus  at  the  top  of  the  clay  3o/3b 
depends  both  on  b  and  t.  At  any  other  height  a,  or  when  the 
conditions  of  simple  vertical  flow  of  pore-water  no  longer  obtain  or  H 
varies  with  time,  the  term  3a/3b  can  be  expected  to  depend  upon 
(a,b,t).  We  shall  now  determine  this  relation. 


Integrating  ( B 1 7)  with  rspect  to  a,  and  using  (B5),  it  is  found 


that 


0  =  Pf  / 


a0(b) 


ao(b) 


da 


l+e  -  +  es  /  -nir  +  Pf[H(tH(a  (b),b,t)]  (B22) 

1  eo  s  a  0 


where  the  last  term  represents  the  water  pressure  on  the  surface 
a  =  a^(b)  of  the  mound,  which  equals  the  total  stress  o  in  (B22)  when 
a  =  aQ(b) . 

From  (B5): 


5(a,b,t) 


a 


/ 

0 


l+e(a,b,t)  , 

i>en(a;br da 


(B23) 


since  5(0, b,t)  =  0.  Substituting  (B23)  into  (B22)  and  rearranging,  it 
is  found  that* 


a  =  p  H  -  pf 


i 3  /  l+e  N  ^  da 

J  TZr~  da  +  (P  ~  Pf)  f  TI —  • 
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(B24) 


It  follows  from  (B24)  that 
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The  first  term  represents  the  weight  of  water  if  it  completely 
filled  the  cylinder  from  the  base  to  the  water  surface.  The  second  term 
is  the  weight  of  water  if  it  filled  the  cylinder  from  the  base  to  the 
height  5(a,b,t)  at  which  a  is  acting.  The  last  term  is  the  buoyant 
weight  of  solids  above  the  plane  on  which  a  acts. 
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const.  to 


which  reduces  in  the  case  = 

(P  ~PJ  dart  p  a  . 

3cj  _  sf  0  f  r  9e 

9b  "  l+eQ  db  "  l+e0  1  9b  a  * 

the  first  term  arising  from  the  initial  non-uniform  surface  profile  of 
the  mound,  the  second  being  a  time-varying  correction  to  this  resulting 
from  the  settlement. 

5.  Some  Typical  Boundary  Conditions 

If  a  mass  of  saturated  clay  of  more-or-less  uniform  void  ratio  e^ 
is  deposited  in  a  short  period  of  time  beneath  water  and  rests  on  (say) 
an  impervious  base,  the  void  ratio  distribution  in  the  mound  at 
subsequent  times  will,  on  the  basis  of  the  assumptions  that  have  been 
made  herein,  be  governed  by  (B19)  with  (B26)  and  subject  to: 

(a)  the  initial  condition 


e(a,b,0)  =  e^  (constant); 

(b)  the  mound  surface  boundary  condition 


e(aQ(b),b,t)  =  eg  (constant); 

(c)  the  impervious  base  condition 


w  (0  ,b,t)  =  0 
a 


v  (0  ,b ,  t )  =  0 
a 


which  result  in  the  following  condition  on  the  void  ratio: 

9e  (ps“Pf)  de 


APPENDIX  C 


In  this  Appendix  the  equations  governing  the  consolidation  of  a 
submarine  clay  mound  under  conditions  of  axial  symmetry  are  derived  in 
Lagrange  coordinates.  The  work  described  in  Appendix  B  was  confined  to 
plane  flow  of  pore-water  and  would  approximate  the  conditions  encountered 
in  long  mounds  of  uniform  cross-section  away  from  the  ends.  The  present 
treatment  extends,  therefore,  the  study  to  mounds  of  compact  form  which 
can  be  approximated  by  bodies  within  a  surface  of  revolution  about  a 
vertical  axis. 

The  physical  assumptions  and  notation  adopted  in  Appendix  B  wi 1 1  be 
retained  here  and  the  development  below  should  be  read  in  coniunction 
with  that  work  upon  which  it  relies  heavily. 

I .  The  Physical  Equations 

The  initial  coordinates  of  a  point  of  the  soil  skeleton  are  (a,b) 
where  £  is  taken  as  vertically  upward  against  gravity  and  b  is  now  a 
radial  distance  measured  horizontally  from  the  axis  of  symmetry  of  the 
mound . 

The  continuity  equation  of  the  solid  phase  retains  the  form 


H  =  1+e 
3a  1+e 

c 

but  that  of  the  pore  water  now  becomes 


^n(va“wa)l  +  fb  [nVb  H]  +  k  [n  If] 


b  3a 


The  flow  of  the  pore  water  is  still  governed  by  the  pair  of 
equations 


(C3) 

(C4) 
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The  spatial  gradients  of  the  excess  pore  water  pressure  (u)  are 
related  to  those  of  the  pore  water  pressure  (p)  through  the  equations 

1“  l  =  Ip  (  °.  \  n  < 


3a  \  l +e  /  3a  \ 1 +e 

3u  _  _3p 
3b  "  3b 


which  allow  (05)  to  be  written  in  terms  of  the  por^  water  pressure  and 
the  void  rat i o  (e ) . 

When  the  equation  of  vertical  equilibrium,  namely 


+  [nPf  +  ( »-n>P8]  H  *  n  . 


is  made  use  of,  together  with  the  definition  of  the  vertical  effective 


stress 


o'  =  o  -  p  , 

equation  (C5)  can  be  written  in  terms  of  the  void  ratio  (e)  above  as 


v  s  '  de  I  1+e  I  3a  1+e  3t 

o 

which  is  the  axisymmetric  Form  of  equation  (BI9)  and  the  governing 
equation  we  seek. 


In  all  other  respects  the  results  arrived  at  for  the  case  of  plane 
pore  water  flow,  and  documented  in  Appendix  B,  hold  mut at i s  mut and i s  for 
this  case  also. 


APPENDIX  D 


In  this  Appendix  we  formulate  the  basic  procedure  for  the  method  of 
lines.  This  is  the  method  we  suggest  for  the  solution  of  the  governing 
equations  presented  in  the  previous  appendices. 

We  consider  a  governing  equation  of  the  form  given  by  equation 
(B19).  In  a  general  form  this  equation  can  be  written  as 

A(e)  ~  +  B(e)(|?-)  +  C(e)  +  D(e)(||)  +  E(e)  —  =  —  (Dl) 

~  2  ^da'  .,2  ^db'  3a  dt 

da  db 

The  method  of  lines  consists  of  a  procedure  in  which  all  but  one  of 
the  independent  variables  are  discretized.  This  yields  a  coupled  system 
of  ordinary  differential  equations.  For  the  equation  given  by  (Dl)  we 
choose  to  discretize  a  into  m  equidistant  mesh  points,  having 
coordinates  (a.:  i  =  l,2,...,m)  and  b  into  n  equidistant  mesh  points, 
having  coordinates  (b . :  i  =  l,2,...,n). 

The  partial  derivatives  are  approximated  at  the  mesh  points  by 

(Zjl)  =  lltlL-Cltld  (D2) 

V  a/i ,  j  2(Aa) 

(-.2  x  e.  .  .-2e .  ,+e.  . 

-ll  -  —  (03) 

3aVi,.i  (Aa)Z 

with  similar  expressions  at  (i)  and  (j-1),  ( i )  and  (i+l)  for  the 
derivatives  with  respect  to  b. 

The  discretized  system  of  governing  equations  now  becomes 
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where  i  =  1,2,... ,m  and  j  =  l,2,...,n.  Thus  we  have  a  system  of  raxn 
ordinary  differential  equations. 

1 .  Method  of  Solution 

The  system  of  ordinary  differential  equations  (D4)  is  solved  by  an 
Adams-Bashforth  predictor  which  is  used  in  conjunction  with  an 
Adams-Moulton  corrector  (26). 

For  the  purpose  of  developing  the  formulation  of  the  method  we 
concentrate  on  a  single  equation  of  the  system  given  by  equation  (D4)  . 

We  aim  to  obtain  the  solution  to 

||=f(t,e(t))  (D5 

in  the  closed  interval  [0,tr.  1.  where  f(t,e(t))  denotes  the  function 

representing  the  left  hand  side  of  the  chosen  equation  from  the  system 

(D4).  To  accomplish  our  objective  we  divide  the  closed  interval 

into  a  set  of  time  steps  (tA,t ,,..., t  )  such  that  (t„  =  0)  and 

0  l  p  0 

(t  =  tc.  ).  We  obtain  the  solution  at  these  time  steps.  The  basic 
p  tin 

computational  task  relates  to  the  advancement  of  the  numerical  solution 
to  (tk+j)  after  having  computed  the  value  of  (e)  at  time  (t^). 

We  note  that 


S<tIr+l  >  =  e(tv)  +  /  af  dt 


or  from  equation  (D5) 


e(tW  =  e(tk)  +  J 

^1. 


f(t,e(t))dt 


The  Adams  method  approximates  the  solution  by  replacing  f(t,e(t)) 
with  a  polynomial  of  order  (i)  which  interpolates  the  computed 

derivatives  at  the  preceding  points  from  the  set  of  discrete  time  values 
and  then  by  integrating  the  polynomial.  Thus  the  polynomial  P.  (t) 

j  K 

must  satisfy  the  following  condition 


P„  ,  (t. )  -  f(t. ,e(t. ));  i  *  0,1,2 . k 

£,k  ill 


The  expression 


ek+l  “  e(V  +  l  k+1  P4,k(t)dt 

fck 

defines  the  order  predictor  of  the  solution  e(t  )  at  time 

(tk+P 

By  choosing  a  constant  time  step 


tk  =  tk.j  +  h  •  k  “  1 . 2  ,  .  .  .  ,  p 


the  backward  formulation  of  the  predictor  is  given  by 


ek+i  3  e(tk>  +  h  J.  Vi?l~  fk 
1*1 


(DIO) 


where 


fk  "  f(tk’  e(tk)) 


1  f  k+l  (t-tk)(t-tk-l)--(t-tk-l-i) 

Yi=T!hJ  - 7^1 - 

ck  h 

and  the  backward  difference  operator  7l  is  defined  as 

Vlf,  =  V( )  =*  Vl-If,  -  7l_If  , 
k  k  k  k-1 


where  we  note  that 


*  fk  3  fk 


7  f,  *  7  f  =  f.  -  f.  , 
k  k  k  k-1 


(Dll) 


(Dl  2) 


(D13) 


(Dl  4) 


(Dl 5a) 


(Dl5b) 


For  purposes  of  efficiency  with  respect  to  the  time  to  solve  the 
system  of  differential  equations  the  Adams  method  uses  a  variable  time 
step  (h^)  and  a  variable  order  of  the  predictor.  The  backward  difference 
formulation  is  replaced  by  a  modified  divided  difference  formulation. 
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as  computed  above  we  evalute  its 


G  iven 
der ivat ive 
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with 


prediction 
respect  to  time 


k+ 1 


h.k" 


e 

k+1 


k+l 
)  »  f 


by 
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(D16) 


The  corrector  to  solve  for  e(t,  ,)  is  formulated  as 

k+ 1 

'"k,i>  -«V  <D17) 

h 

The  polynomial  P^+^  ^Ct)  interpolates  the  same  data  as  P^ 
with  an  additional  point,  namely 


Pl+1 ,k(tk+l)  ~  fk+l 


(D18) 


It  now  remains  only  to  evaluate  f(t^+p  e(t^+j)  for  the  next  time 

step . 

The  algorithm  described  (26)  uses  low  orders  of  the  predictor 
polynomials  (P).  This  maximizes  the  stability  properties.  The  order  of 
the  predictor  polynomial  is  accomplished  separately  from  the  time  step 
selection. 


The  time  step  is  determined  as  the  largest  value  for  which  the  local 
error  meets  a  given  input  tolerance.  The  control  of  the  local  error 
assumes  that  the  local  solution  is  uniformly  approximated  over  the  closed 
interval  [t^,  t^+j].  The  local  error  is  controlled  per  unit  time  step; 
i.e.  relative  to  the  size  of  the  step.  Using  this  procedure  the  solution 
overshoots  the  last  time  ( t . n ) .  Since  the  approximation  is  uniform 
over  the  whole  interval,  however,  an  interpolation  provides  the  solution 
at  the  required  value  of  time. 

The  procedure  sets  all  the  parameters  for  continuation  of  the 
predictor-corrector  process.  Thus  a  computation  at  a  predetermined 
sequence  of  time  is  provided  by  this  algorithm. 


